#TEST#
#' Detect cloud, cloud shadow and snow using the Fmask method
#'
#' The Fmask (Function of mask) is an automated cloud, cloud shadow and snow masking detect ion method for Landsat TM/ETM+ images developed by Zhu and Woodcock (2012). A Matlab code for the Fmask is maintained by Zhu (http://code.google.com/p/fmask/). This is translation of the Matlab code to R. 
#'
#' @name FmaskforR-package
#' @aliases FmaskforR
#' @docType package
#' @title Automated cloud, cloud shadow and snow masking detection method for Landsat TM/ETM+ images.
#' @author Joseph Henry \email{joseph.henry@@sydney.edu.au} and Willem Vervoort \email{willem.vervoort@@sydney.edu.au}
#' @references Zhu, Z. and Woodcock, C.E. 2012. Object-based cloud and cloud shadow detection in Landsat imagery. Remote Sensing of Environment 118, 83-94. ' \url{http://code.google.com/p/fmask/}
#' @keywords package, Fmask, Landsat, Cloud, Snow, Shadow

### CONDITION OF USE ###
# Joseph and Willem reserve the right to use this code or its derivatives in the creation of an R package.

### ALPHA VERSION ###

	## Packages (dependicies of the package)
	 require (stringr)
	 require (raster)
	 require(RPyGeo)
	 

## ArcGIS setup
fill_func <- function (in_surface_raster, out_surface_raster, z_limit) {
	rpygeo.geoprocessor (ARCGIS.FILLSA,
	c (in_surface_raster, out_surface_raster, z_limit),
	env = env.cur,
	 working.directory = reproj_d)
	}
## Configure python
  env.cur <- rpygeo.build.env(
  workspace = "C:/",# ARC related
  python.path = "C:/Python25",python.command = "python.exe" #python related
  ,overwriteoutput = 1 #if 1, overwrite
  )
	# Working directory
   reproj_d <- paste("C:/",sep="")


shad_cell_func <- function (input.a, input.b){
#' Function used in Fmask.
	shadow_cells <- cellFromRowCol (input.a, input.b[,1], input.b[,2])
	return (shadow_cells)
}

FmaskForR <- function (INPUT_BANDS, INPUT_BANDS_DN, run.prob.over.water, L7.SUN.AZIMUTH.TAN.RADIANS, L7.SUN.AZIMUTH, L7.SUN.ELEVATION.TAN.RADIANS,OUTPUT) {
#' Produce Landsat reflectance and temperature rasters.
#' @param INPUT_BANDS file paths of reflectance and temperature data (bands 1 to 7) (character)
#' @param INPUT_BANDS_DN file paths of digital number data (bands 1 to 3) (character)
#' @param run.prob.over.water choice,  "yes" or  "no", to run Fmask over water. JH didn't use this because no large water bodies in catchment.
#' @param OUTPUT file path of output file (character).
#' @param L7.SUN.AZIMUTH.TAN.RADIANS sun azimuth in radians derived from metadata file that comes with tile (numeric)
#' @param L7.SUN.AZIMUTH sun azimuth in degrees derived from metadata file that comes with tile (numeric)
#' @param L7.SUN.ELEVATION.TAN.RADIANS  sun elevation derived from metadata file that comes with tile (numeric)
 ## Load
#################Constants set for testing 
  # INPUT_BANDS <- DATA_TR_OUT
	# INPUT_BANDS_DN <- DATA_PATH_RAW_EXCL
	# OUTPUT <- "C:/fmask.tif"
#################
	# Rasters
	 # Reflectance and temperature
    Band_1 <- raster (INPUT_BANDS [1]) # Band 1	
	  Band_2 <- raster (INPUT_BANDS [2]) # Band 2		
	  Band_3 <- raster (INPUT_BANDS [3]) # Band 3		
	  Band_4 <- raster (INPUT_BANDS [4]) # Band 4		
	  Band_5 <- raster (INPUT_BANDS [5]) # Band 5		
	  Band_6 <- raster (INPUT_BANDS [6]) # Band 6 in Celcius		
	  Band_7 <- raster (INPUT_BANDS [7]) # Band 7	
	 # Digital numbers
		Band_1_DN <- raster (INPUT_BANDS_DN [1])[] # Band 1	
		Band_2_DN <- raster (INPUT_BANDS_DN [2])[] # Band 2
		Band_3_DN <- raster (INPUT_BANDS_DN [3])[] # Band 3		
	# Vectors
	 Band_1_vec <- Band_1[]
	 Band_2_vec <- Band_2[]
	 Band_3_vec <- Band_3[]
	 Band_4_vec <- Band_4[]
	 Band_5_vec <- Band_5[]
	 Band_6_vec <- Band_6[]
	 Band_7_vec <- Band_7[]
	rm (Band_2); rm (Band_3); rm (Band_4); rm (Band_5); rm (Band_6); rm (Band_7); # Remove some of the rasters

	## Clean
	 # Identify pixels that are saturated in the visible L7.PASS.CHOICE. Zhu has value of 255
	  B1Satu <- which (Band_1_DN == L7.FILLVALUE)
		B2Satu <- which (Band_2_DN == L7.FILLVALUE)
		B3Satu <- which (Band_3_DN == L7.FILLVALUE)
		B_Satu_unique_all <- sort (unique (c (B1Satu, B2Satu, B3Satu)))
		B_Satu_table <- table( as.numeric (c (B1Satu, B2Satu, B3Satu)))
		B_Satu_rows_all <- as.numeric (names (B_Satu_table[which (B_Satu_table == 3)]))
	 # Missing values
	  id_missing <- unique(c(
		which(Band_1_vec_raw_L7.FILLVALUE == min (Band_1_vec_raw_L7.FILLVALUE)),
		which(Band_2_vec_raw_L7.FILLVALUE == min (Band_2_vec_raw_L7.FILLVALUE)),
		which(L7.VECTOR.BAND.3.RAW_L7.FILLVALUE == min (L7.VECTOR.BAND.3.RAW_L7.FILLVALUE)),
		which(L7.VECTOR.BAND.4.RAW_L7.FILLVALUE == min (L7.VECTOR.BAND.4.RAW_L7.FILLVALUE)),
		which(Band_5_vec_raw_L7.FILLVALUE == min (Band_5_vec_raw_L7.FILLVALUE)),
		which(Band_6_vec_raw_L7.FILLVALUE == min (Band_6_vec_raw_L7.FILLVALUE)),
		which(Band_7_vec_raw_L7.FILLVALUE == min (Band_7_vec_raw_L7.FILLVALUE))))
	 # Convert 0 to NA.NUMERIC
		Band_1_vec_na <- (Band_1_vec)
		Band_2_vec_na <- (Band_2_vec)
		Band_3_vec_na <- (Band_3_vec)
		Band_4_vec_na <- (Band_4_vec)
		Band_5_vec_na <- (Band_5_vec)
		Band_6_vec_na <- (Band_6_vec)
		Band_7_vec_na <- (Band_7_vec)
		if (length (id_missing) > 0){ # Set NAs
			Band_1_vec_na[id_missing] <- NA
			Band_2_vec_na[id_missing] <- NA
			Band_3_vec_na[id_missing] <- NA
			Band_4_vec_na[id_missing] <- NA
			Band_5_vec_na[id_missing] <- NA
			Band_6_vec_na[id_missing] <- NA
			Band_7_vec_na[id_missing] <- NA
		}	
   # Remove some of the rasters
	  rm (Band_1_DN); rm (Band_2_DN); rm (Band_3_DN);	

  ## Do		
	 # If all values are not NA
	  if (length (id_missing) != length (Band_7_vec)){
		 satu_Bv.data <- (Band_1_vec_na + Band_2_vec_na + Band_3_vec_na)[B_Satu_unique_all]
		 satu_Bv <- satu_Bv.data >= 1
		 Satu_rows <- B_Satu_unique_all[which (satu_Bv == TRUE)]
    # Potential cloud layer -- pass one
  	 BT <- Band_6_vec_na # Renamed (same as matlab name)
			## Eq 1: Basic cloud test
			 NDVI <- (Band_4_vec_na - Band_3_vec_na) / (Band_4_vec_na + Band_3_vec_na) # Normalised Difference Vegetation Index
			 NDSI <- (Band_2_vec_na - Band_5_vec_na) / (Band_2_vec_na + Band_5_vec_na) # Normalised Difference Snow Index						 
			 Basic_Test <- ifelse (Band_7_vec_na > 300 & BT < 2700 & NDSI < 0.8 & NDVI < 0.8, TRUE, FALSE) # Run test
			## Eq 2: Whiteness test
			 MeanVis <- (Band_1_vec_na + Band_2_vec_na + Band_3_vec_na) / 3
			 B1.MeanVis.diff <- abs (Band_1_vec_na - MeanVis)
			 B2.MeanVis.diff <- abs (Band_2_vec_na - MeanVis)
			 B3.MeanVis.diff <- abs (Band_3_vec_na - MeanVis)
			 Bsum.MeanVis.diff <- (B1.MeanVis.diff + B2.MeanVis.diff + B3.MeanVis.diff) / MeanVis
			 whiteness <- Bsum.MeanVis.diff
			## Adjust whiteness values
			 # Set saturated whiteness to zero
			 if (length (Satu_rows) > 0){
			  whiteness[Satu_rows] <- 0
			 }
			 whiteness[whiteness > 1] <- 1 #C- whiteness goes over 1, but I don't know why. Therefore reset to 1. -C#	
			 Whiteness_Test <- whiteness < 0.7
		  ## Eq 3: Haze test
			 HOT_Test.data <- Band_1_vec_na - 0.5 * Band_3_vec_na - 800 # Test values
			 HOT_Test <- HOT_Test.data > 0 # Run test
			## Eq 4: Band 4 and 5 ratio
       B45_Test <- (Band_4_vec_na / Band_5_vec_na) > 0.75 # Test
      ## Eq 5: Water test
			 Water_Test <- ifelse ((NDVI < 0.01 & Band_4_vec_na < 1100)|(NDVI < 0.1 & NDVI > 0 & Band_4_vec_na < 500), TRUE, FALSE)
      ## Eq 6: Potential cloud ?
			 PCP <- ifelse (Basic_Test == TRUE & Whiteness_Test == TRUE & HOT_Test == TRUE & B45_Test ==TRUE, TRUE, FALSE)
				
			####S Choice switch
				choice_switch_1 <- length (which (PCP == TRUE)) / length (PCP)
			####E
			
				if (choice_switch_1 >= 0.99) { # Too many clouds, go straight to end. Set entire tile to cloud.
				 Fmask <- rep (1, length (Band_1_vec_na)) # Set vector
				 Fmask_ras <- setValues (Band_1, Fmask) # Turn vector into raster
				 writeRaster(Fmask_ras,fmask.dir) # Save raster
				}
					
				if (choice_switch_1 < 0.99){ # Not fully cloud. Potential cloud layer -- pass two 
				# Temperature probability for water
				 if (run.prob.over.water == "yes"){ # Are there any bodies of water to run this over? 
          ## Eq 7: Temperature test over water
					 Clearsky_Water <- (Water_Test == TRUE & Band_7_vec_na <= 300) #C- Inconsistency between matlab and paper? -C#			 
					## Eq 8
					 T_water <- quantile (as.numeric (BT[Clearsky_Water == TRUE]), 0.825, na.rm = TRUE)
					## Eq 9
					 wTemperature_Prob <- (T_water - BT) / 400
					## Eq 10: Brightness test over water
					 t_bright <- 1100
					 Brightness_Prob <- Band_5_vec_na / t_bright
					 Brightness_Prob[Brightness_Prob > 1] <- 1
					## Eq 11: Final probability mask over water
					 if (is.na (T_water) == FALSE){
					  wCloud_Prob <- wTemperature_Prob * Brightness_Prob * 100
						} else { #C- If no clear-sky water was identified -C#
						wCloud_Prob <- 1 / (Brightness_Prob) 
					 }
				 }
				# Temperature probability for land
				 ## Eq 12
					T_buffer <- 4 * 100
					Clearsky_Land <- ifelse (PCP == FALSE & Water_Test == FALSE, TRUE, FALSE)
				 ## Eq 13	
				  T_low <- quantile (as.numeric (na.omit (BT[Clearsky_Land == TRUE])), 0.175) # Low percentile
					T_high <- quantile (as.numeric (na.omit (BT[Clearsky_Land == TRUE])), 0.825) # High percentile
					# Other stuff
						t_tempL <- T_low - T_buffer
						t_tempH <- T_high + T_buffer
						Temp_l <- t_tempH - t_tempL
				 ## Eq 14
				  lTemperature_Prob <- (t_tempH - BT) / Temp_l
					lTemperature_Prob[lTemperature_Prob < 0] <- 0 #C- Probability can be > 1 -C#
					
				####S Choice switch
					choice_switch_2 <- length (which (Clearsky_Land == TRUE)) / length (Clearsky_Land)
				####E
				}
				
				#C- Fmask will use the clear clear-sky pixels from both land and water for calculating temperature probability -C#
				if (choice_switch_1 < 0.99){
					if (choice_switch_2 < 0.001){ 	
					# NDSI
						modified_NDSI <- NDSI
						modified_NDSI_tf_1 <- NDSI[B3Satu]<0
						modified_NDSI_tf_2 <- B3Satu[which(modified_NDSI_tf_1==TRUE)]
						modified_NDSI[modified_NDSI_tf_2] <- 0								
					## Eq 20: potential snow layer
					 PSL <- ifelse(modified_NDSI>0.15&BT<380&Band_4_vec_na>1100&Band_2_vec_na>1000,TRUE,FALSE)
					# Vector with values 4) cloud pixels, 2) shadow pixels, 3) snow pixels
 					 Fmask <- rep(0,length(Band_1[]))
					 Fmask[PSL == TRUE] <- 4 	# Snow
					 Fmask[id_missing] <- NA # NA post 2003
  				 Fmask_ras <- setValues(Band_1,Fmask) # Create raster
					 writeRaster(Fmask_ras,fmask.dir) # Save raster
					}
				}
			
			#C- Fmask will use the clear clear-sky pixels land only for calculating temperature probability -C#
				if (choice_switch_1 < 0.99){				
					if (choice_switch_2 >= 0.001){			
						## Variability probability
						 # Modify vales for saturated pixels
						  # NDVI
							 modified_NDVI <- NDVI								
							 modified_NDVI_tf_1 <- NDVI[B2Satu] > 0
							 modified_NDVI_tf_2 <- B2Satu[which (modified_NDVI_tf_1 == TRUE)]
							 modified_NDVI[modified_NDVI_tf_2] <- 0
							# NDSI
							 modified_NDSI <- NDSI
							 modified_NDSI_tf_1 <- NDSI[B3Satu] < 0
							 modified_NDSI_tf_2 <- B3Satu[which (modified_NDSI_tf_1==TRUE)]
							 modified_NDSI[modified_NDSI_tf_2] <- 0								
							## Eq 15
#START###################################################
	# COMMENT
	# JH (25/4/12): Change from the Matlab or paper.					
								#C- I removed whiteness because it can equal 1. Hence, Variability prob would be 0. Therefore, subsequent probabilities would be 0 -C#
								Variability_Prob <- 1 #-max(abs(modified_NDVI),abs(modified_NDSI),na.rm=T)#,whiteness,na.rm=T)
								#Variability_Prob <- 1-max(abs(modified_NDVI))#,abs(modified_NDSI))
#END#####################################################							
						## Cloud probability for land
						 ## Eq 16
						  lCloud_Prob <- (lTemperature_Prob * Variability_Prob) * 100
						 ## Eq 17
							 cldprob <- 22.5 # Default cloud probability threshold for cloud detection
               Land_threshold <- as.numeric (quantile (lCloud_Prob[Clearsky_Land == TRUE], 0.825, na.rm = TRUE) + cldprob) 	# Treshold for defining cloud cover
					  ## Potential cloud layer
						 ## Eq 18
						  if (run.prob.over.water=="yes"){							
							 clr_max <- 50 # Treshold
							 PCL_t_1 <- PCP == TRUE & Water_Test == TRUE & wCloud_Prob > clr_max # Cloud over water
							}
 							PCL_t_2 <- PCP== TRUE & Water_Test == FALSE & lCloud_Prob > Land_threshold # Cloud over land
							PCL_t_3 <- lCloud_Prob > 99 & Water_Test == FALSE # High probability of cloud over land
							PCL_t_4 <- BT < (T_low - 3500) # Really cold cloud
						## Classify each pixel
							PCL_raw_class <- rep (NA, length (PCL_t_2)) # Blank vector
							# Add values to vector
							 if (run.prob.over.water=="yes"){ PCL_raw_class[PCL_t_1] <- 1 }				  							
							 PCL_raw_class[PCL_t_2] <- 2
							 PCL_raw_class[PCL_t_3] <- 3
							 PCL_raw_class[PCL_t_4] <- 4
							 PCL_raw_class[is.na(PCL_raw_class)] <- 0
							# Make a 0,1 mask		
								PCL <- PCL_raw_class
								PCL[PCL!=0] <- 1
						## Release memory
							if  (run.prob.over.water == "yes"){ rm (PCL_t_2)}
							rm (PCL_t_2); rm (PCL_t_3); rm (PCL_t_4);
						## Cloud work
						 # Clumping (group unique clouds)
							PCL_ras <- setValues(Band_1,PCL)
							PCL_clump <- clump (PCL_ras)
							PCL_clump_vec <- PCL_clump[]
						## Add 8 pixels around each cloud pixel in a clump
						 Cloud_pix_row <- which (PCL == 1) # Cloud pixels
						# Adjacency (outside the edge of cloud)
							PCL_adj_1 <- adjacency (PCL_ras, Cloud_pix_row, 1 : ncell(PCL_clump), directions = 8)
							PCL_adj_v1 <- vector (length = ncell (PCL_clump))			
							PCL_adj_v1[PCL_adj_1[,2]] <- 2
							PCL_adj_v1[Cloud_pix_row] <- 1		
						# Adjacency (at the edge of a cloud)
							Outer_edge_cloud_pix <- which (PCL_adj_v1 == 2)
							PCL_adj_2 <- adjacency (PCL_ras, fromCells = Outer_edge_cloud_pix,toCells = 1:ncell(PCL_clump), directions = 8)
							PCL_adj_v2 <- vector(length = ncell (PCL_clump))			
							Inner_edge <- which (PCL_adj_2[,2] %in% Cloud_pix_row)
							PCL_adj_v2[PCL_adj_2[,2]] <- 3
							PCL_adj_v2[PCL_adj_v2 != 3] <- 0
							PCL_adj_v2_2 <- rep (NA, length (PCL_adj_v2))
							PCL_adj_v2_2[PCL_adj_v1==1&PCL_adj_v2==3] <- 1
							PCL_adj_v2_2[is.na(PCL_adj_v2_2)==TRUE] <- 0
							Cloud_edge_pix <- PCL_adj_v2_2
						# Re-clump without edges
							PCL_edge <- PCL
							PCL_edge[which (Cloud_edge_pix == 1)] <- 0
							# Can't really clump if only 1 pixel
								if (length(PCL_edge[PCL_edge==1]) > 1){
								 PCL_ras_edge <- setValues (Band_1,PCL_edge)
								 PCL_clump_edge <- clump (PCL_ras_edge)
								 PCL_clump_vec_edge <- PCL_clump_edge[]
								 } else {
								 PCL_clump_edge <- PCL_ras_edge
								 PCL_clump_vec_edge <- PCL_edge
								}
					## EQ 22/23 Cloud base 
 					 BT_split <- split (BT, PCL_clump_vec_edge) # Split BT values
					 # Length of clouds
					  PCL_clump_vec_edge.table <- table (PCL_clump_vec_edge)
					  PCL_clump_vec_edge.table.df <- data.frame (VAL = names (PCL_clump_vec_edge.table), COU = as.numeric (PCL_clump_vec_edge.table))
						Tcloud_base <- rep (NA, nrow (PCL_clump_vec_edge.table.df)) # Blank vector
					# Cloud base BT
					 for(a in 1 : length(Tcloud_base)){
					  BT_split_iter <- BT_split[[a]]
						Cloud_id <- PCL_clump_vec_edge.table.df[a,1]
						R <- sqrt (PCL_clump_vec_edge.table.df[a,2] / (2 * pi))	
						 if (R >= 8){ #C- Between 0 and 1 in R -C#
						  quant.val <- (((R - 8)^2)) / R^2
							Tcloud_base[a] <- as.numeric (quantile (BT_split_iter,quant.val, na.rm = TRUE))
						 }
						 if (R < 8){
						  Tcloud_base[a] <- min (BT_split_iter)
							Tcloud_base[a] <- ifelse (Tcloud_base[a] < 0, NA, Tcloud_base[a]) #c- read through more -c#
						 }
					 }
					## EQ 23: Tcloud_object > base == base
					 BT_split_mod <- BT_split
					 for(a in 1:length (Tcloud_base)){
					  in_1 <- BT_split[[a]]
					  in_2 <- Tcloud_base[a]
						in_1[in_1 > in_2] <- in_2
						BT_split_mod[[a]] <- in_1
					 }
					 BT_mod <- unsplit (BT_split_mod, PCL_clump_vec_edge) # Modified BT
				  ## Eq 20: potential snow layer
					 PSL <- ifelse (modified_NDSI > 0.15 & BT < 380 & Band_4_vec_na > 1100 & Band_2_vec_na > 1000, TRUE, FALSE)
			    ## Object-based cloud and cloud shadow match
				   revised_ptm <- 100 * (length (which (PCL == 1)) / length (PCL))
					 if (revised_ptm >= 90){ # Too many clouds.  Set entire layer to cloud.
					  print("No shadow match because too much cloud")
						Fmask <- rep (1, length (Band_1_vec_na)) # Set vector values
						Fmask_ras <- setValues (Band_1, Fmask)
						writeRaster (Fmask_ras, fmask.dir)
				   }
					if (revised_ptm < 90){ # Not too many clouds
#START###################################################
	# COMMENT
	# JH (25/4/12): Can't get TauDEM to run. Script below is for ArcGIS. File saved to work directory.
	# JH (25/4/12): Unable to do using TauDEM on my computer so am instead using of ArcGIS. Need to also install http://www.mcs.anl.gov/research/projects/mpich2/downloads/index.php?s=downloads.
	# JH (25/4/12): I don't know how to do this in R. I would like to because it would avoid the need for Python and ArcGIS installation.
						Band_4_vec_ref <- setValues (Band_1, as.numeric (Band_4_vec_na))
						writeRaster(Band_4_vec_ref,  "BAND_4_TEMP.tif", overwrite = TRUE)
						fill.mv <- max (Band_4_vec_na, na.rm = TRUE)
						fill_func (paste (getwd(), "BAND_4_TEMP.tif", sep = "/"), paste (getwd(), "BAND_4_TEMP_FILL.tif", sep = "/"), NULL)
						Band_4_filled_ras <- raster ("BAND_4_TEMP_FILL.tif")
						Band_4_filled_vec_na <- Band_4_filled_ras[]
						rm (Band_4_filled_ras)
           # Repeat for Band 5
						Band_5_vec_ref <- setValues (Band_1, as.numeric (Band_5_vec_na))
						writeRaster(Band_5_vec_ref,  "BAND_5_TEMP.tif", overwrite = TRUE)
						fill.mv <- max (Band_5_vec_na, na.rm = TRUE)
						fill_func (paste (getwd(), "BAND_5_TEMP.tif", sep = "/"), paste (getwd(), "BAND_5_TEMP_FILL.tif", sep = "/"), NULL)
						Band_5_filled_ras <- raster ("BAND_5_TEMP_FILL.tif")
						Band_5_filled_vec_na <- Band_5_filled_ras[]
						rm (Band_5_filled_ras)
#END#####################################################	
					## Eq 19
					 shadow_prob <- apply (cbind (Band_4_filled_vec_na - Band_4_vec_na,Band_5_filled_vec_na - Band_5_vec_na), 1, function(x) min (x, na.rm = TRUE))
					 PCSL <- shadow_prob > 200	
#START###################################################
	# COMMENT
	# JH (25/4/12): This the iteration part that really slows the process down. Parallel computing here would be good.	
					## Constants
					 t_obj <- Tcloud_base
					 Tsimilar <- 0.3
					 Tbuffer <- 0.98 # Treshold for matching buffering					
					 num_cldoj <- 25 # Minimum matched cloud object (pixels)
					 num_pix <- 8 # Number of inward pixels (240 M) for cloud base temperature
					 rate_elapse <- 6.5 # Lapse rates
					 rate_dlapse <- 9.8 # Lapse rates
					 Max_cl_height <- 12000 	# Max cloud base height (m)
					 Min_cl_height <- 200 # Min cloud base height (m)				
					## Cloud heights
					 Min_cl_height_1 <- unlist (lapply (t_obj, function (x) max (Min_cl_height, 10 * (t_tempL - 400 - x) / rate_dlapse, na.rm = TRUE)))
					 Min_cl_height_1 <- ifelse (Min_cl_height_1 < 0, NA, Min_cl_height_1)
					 Max_cl_height_1 <- unlist (lapply (t_obj, function (x) min (Max_cl_height, 10 * (t_tempH - 400 - x), na.rm = TRUE)))
					 Max_cl_height_1 <- ifelse (Max_cl_height_1 < 0, NA, Max_cl_height_1)
				 ## Object-based cloud and cloud shadow match
					BT_mod_split <- split (BT_mod, PCL_clump_vec_edge) #C- Iteration for each cloud -C#
					PCL_clump_vec_edge.mem <- data.frame (ID = 1:length (PCL_clump_vec_edge), PCL_clump_vec_edge = PCL_clump_vec_edge) #C- made due to memory problems -C#					
					shadow_cells_list <- list()	
					h <- list()
					YX_mod_out_list <- list()
					## Each cloud
					 for(a in 1:length (Tcloud_base)){
					  if (is.na (Min_cl_height_1[a]) == FALSE & is.na (Max_cl_height_1[a]) == FALSE){
						 shadow_cells_list_out <- list() # List
						 sub_size <- xres (Band_1) 	# Resolution (assumes xres==yres)
						 #i_step <- 2* # Move 2 pixels at a time
						 # Calculate cloud DEM with initial base height
						  base_h <- seq (Min_cl_height_1[a], Max_cl_height_1[a], 10)
							t_obj_a <- t_obj[a]
						 #base_h <- Min_cl_height_1[a]
						 temp_obj <- BT_mod_split[[a]]
						 ## For each height of cloud
						 for(b in 1:length (base_h)){
							h <- 10 * (t_obj_a-temp_obj) / rate_elapse + base_h[b] # Each height in the cloud
							i_xy <- h / (sub_size * tan (L7.SUN.ELEVATION.TAN.RADIANS)) 	# Shadow moved distance (pixel)
							# Movement (pixels)
								X_move <- round ((i_xy * cos (L7.SUN.AZIMUTH.TAN.RADIANS)) / 30)
								Y_move <- round ((i_xy * sin (L7.SUN.AZIMUTH.TAN.RADIANS)) / 30)
							# Indentify cloud cells
							 cloud.cell.id <- as.numeric (na.omit (
							 PCL_clump_vec_edge.mem[PCL_clump_vec_edge.mem[,2] == a, 1]))
								#which(PCL_clump_vec_edge==a)	
							 YX_ori <- rowColFromCell (PCL_clump_edge, cloud.cell.id)	
							# Identify shadow rowCol positions
							 YX_mod <- YX_ori
							 if (L7.SUN.AZIMUTH < 180){
							 #if (L7.SUN.AZIMUTH.TAN.RADIANS < 180){
							  YX_mod[,1] <- YX_ori[,1] - Y_move
							  YX_mod[,2] <- YX_ori[,2] - X_move
								} else {
								YX_mod[,1] <- YX_ori[,1] + Y_move
								YX_mod[,2] <- YX_ori[,2] + X_move				
							 }
							# Values that are off the map
							 YX_mod[YX_mod[,1] <= 0,] <- NA
							 YX_mod[YX_mod[,2] <= 0,] <- NA
							 YX_mod_out <- na.omit (YX_mod)
							if (b == 1){
							 shadow_cells_list_out_b <- shad_cell_func (PCL_clump_edge, as.data.frame (YX_mod_out, stringsAsFactors = FALSE))
							 } else {
							 shadow_cells_list_out_b <- as.numeric (unique (na.omit (c (shadow_cells_list_out_b, shad_cell_func (PCL_clump_edge, as.data.frame (YX_mod_out, stringsAsFactors = FALSE))))))
							}												
							rm(YX_mod_out); rm(YX_mod); # Free memory
							print(paste(b, length(base_h),sep=" of ")) # Print progress
						 }
  					 } else {
						 shadow_cells_list_out_b <- NA
						 }
						 # Output for each cloud				
						  if (a == 1){
							 shadow_cells_list <- shadow_cells_list_out_b
							 } else {
							 shadow_cells_list <- as.numeric (unique (na.omit (c (shadow_cells_list, shadow_cells_list_out_b))))
							}
							# See where it is at
								print("#####################################")
								print (paste (a, length (Tcloud_base), sep = " of "))
								rm(YX_ori);rm(cloud.cell.id);rm(Y_move);rm(X_move);rm(i_xy);rm(h);rm(t_obj_a);rm(temp_obj);rm(sub_size); rm(shadow_cells_list_out_b); # Free memory
					 }					
					## Cloud raster
					 Unique.cloud.pix <- unique(shadow_cells_list) # Unique pixels only
					 if (length(Unique.cloud.pix)!=0){
						cloud_shadow <- rep (NA, length (PCL_clump_edge[])) # Make vector
						cloud_shadow[unlist (Unique.cloud.pix)] <- 1 # Shadow							
						cloud_shadow_ras_1 <- setValues (Band_1, cloud_shadow) # Make raster
						## Dialte cloud and shadow pixels 3 times. Zhu and Woodcock, 2012 call for buffers by 3 pixels in 8-connected directions -C#
						# 1 pixel out
							shadow_cal_1 <- adjacency (cloud_shadow_ras_1, fromCells = Unique.cloud.pix, toCells = c(1:ncell (cloud_shadow_ras_1)), directions = 8)
							shadow_cal_1_vec <- rep (NA, length (PCL_clump_edge[]))		
							shadow_cal_1_vec[shadow_cal_1[,2]] <- 1
							shadow_cal_1_ras <- setValues (Band_1, shadow_cal_1_vec)				
						# 2 pixels out
							shadow_cal_2_vec_vals <- which (shadow_cal_1_vec == 1)
							shadow_cal_2 <- adjacency (shadow_cal_1_ras, fromCells = shadow_cal_2_vec_vals, toCells = c(1:ncell (shadow_cal_1)), directions = 8)
							shadow_cal_2_vec <- rep (NA, length (PCL_clump_edge[]))		
							shadow_cal_2_vec[shadow_cal_2[,2]] <- 1
							shadow_cal_2_ras <- setValues (Band_1, shadow_cal_2_vec)
						# 3 pixels out
							shadow_cal_3_vec_vals <- which (shadow_cal_2_vec == 1)
							shadow_cal_3 <- adjacency (shadow_cal_2_ras, fromCells = shadow_cal_3_vec_vals, toCells = c(1:ncell (shadow_cal_1)), directions = 8)
							shadow_cal_3_vec <- rep (NA, length (PCL_clump_edge[]))		
							shadow_cal_3_vec[shadow_cal_3[,2]] <- 1
							rm(shadow_cal_2_ras); rm(shadow_cal_3);rm(shadow_cal_3_vec_vals);rm(shadow_cal_2_vec);rm(shadow_cal_2);rm(shadow_cal_2_vec_vals);rm(shadow_cal_1_vec);rm(shadow_cal_1) # Free memory
						} else {
						 shadow_cal_3_vec <- rep (0, length (PCL_clump_edge[]))
					  }

					## Classify output #C- Priority: 1) cloud pixels, 2) shadow pixels, 3) snow pixels
					 Fmask <- rep(0,length(PCL_clump_edge[]))						
					 Fmask[PSL==TRUE] <- 3 # Snow							
					 Fmask[shadow_cal_3_vec[]==1] <- 2 # Shadows								
					 Fmask[PCL==1] <- 1 # Cloud
					 Fmask[id_missing] <- NA # NA post 2003
					 Fmask_ras <- setValues(Band_1,Fmask) # Raster
					 writeRaster(Fmask_ras, OUTPUT, overwrite = TRUE)
						
					## Free memory
					 rm(Fmask_ras); rm(Fmask_ras); rm(id_missing); rm(PCL); rm(shadow_cal_3_vec); rm(PSL); rm(Fmask); rm(PCL_clump_edge);
					 rm(cloud_shadow_ras_1);rm(cloud_shadow);
					 rm(B2.MeanVis.diff); rm(B3.MeanVis.diff); rm(B3.MeanVis.diff); rm(B45_Test); rm(Band_1_vec); rm(Band_1_vec_na); rm(Band_1_vec_raw_L7.FILLVALUE); rm(Band_2_vec); rm(Band_2_vec_na); rm(Band_2_vec_raw_L7.FILLVALUE); rm(Band_3_vec); rm(Band_3_vec_na); rm(L7.VECTOR.BAND.3.RAW_L7.FILLVALUE); rm(Band_4_filled_vec_na); rm(Band_4_vec); rm(Band_4_vec_na); rm(L7.VECTOR.BAND.4.RAW_L7.FILLVALUE); rm(Band_5_filled_vec_na); rm(Band_5_vec); rm(Band_5_vec_na); rm(Band_5_vec_raw_L7.FILLVALUE); rm(Band_6_vec); rm(Band_6_vec_na); rm(Band_6_vec_raw_L7.FILLVALUE); rm(Band_7_vec); rm(Band_7_vec_na); rm(Band_7_vec_raw_L7.FILLVALUE); rm(Basic_Test); rm(Bsum.MeanVis.diff); rm(BT); rm(BT_mod); rm(BT_split); rm(BT_split_mod); rm(Cloud_edge_pix); rm(Cloud_pix_row); rm(HOT_Test); rm(HOT_Test.data); rm(Inner_edge); rm(MeanVis); rm(Outer_edge_cloud_pix); rm(PCL_adj_1); rm(PCL_adj_2); rm(PCL_adj_v1); rm(PCL_adj_v2); rm(PCL_adj_v2_2); rm(PCL_clump); rm(PCL_clump_edge); rm(PCL_clump_vec); rm(PCL_clump_vec_edge); rm(PCL_clump_vec_edge.table); rm(PCL_edge); rm(PCL_ras); rm(PCL_ras_edge); rm(PCP); rm(PSL); rm(satu_Bv); rm(satu_Bv.data); rm(Satu_rows); rm(shadow_prob); rm(Tcloud_base); rm(Water_Test); rm(whiteness); rm(Whiteness_Test); 
				 }
				}
			 }
	    }
}